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Theoretical optical and x-ray spectra of model structures of water and ice are calculated us- 
ing a many-body perturbation theory, Bethe-Salpeter equation (BSE) approach implemented in 
the valence- and core-excitation codes AI2NBSE and OCEAN. These codes use ab initio density- 
functional theory wave functions from a plane-wave, pseudopotential code, quasi-particle self-energy 
corrections, and a BSE treatment of particle-hole interactions. The approach improves upon 
independent-particle methods through the inclusion of a complex, energy-dependent self-energy 
and screened particle-hole interactions to account for inelastic losses and excitonic effects. These 
many-body effects are found to be crucial for quantitative calculations of ice and water spectra. 
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I. INTRODUCTION 

Recently there has been considerable interest and con- 
troversy surrounding the connection betvifcen the local 
structure of water and ice and their observed optical and 
x-ray spectrai^i Part of the difficulty lies in modeling 
the structures of these complex, finite-temperature sys- 
tems. Conventional super-cell methods require large unit 
cells to treat proton-disorder in ice and a configurational 
average adequate to represent a statistical ensemble of 
water structures. Structural probes such as x-ray and 
neutron scattering do not provide an unambiguous in- 
terpretation of local geometry.^ Another part of the dif- 
ficulty lies in theoretical modeling. Despite many at- 
tempts, quantitative theoretical calculations of the opti- 
cal and x-ray spectra of these systems have proved to be 
notoriously difficult due to strong non-local, self-energy, 
and excitonic effects. Thus the various theoretical meth- 
ods that have been employed to date exhibit consider- 
able variation in their results, undermining a definitive 
interpretation!^"— iii^ 

In an effort to address these issues, we present calcu- 
lations of the valence and core excitation spectra of well- 
characterized model ice and water systems based on a 
recently developed approach utilizing the Bethe-Salpeter 
equation (BSE) and Hedin's GW approximation for the 
quasi-particle self-energy (the acronym GW refers to the 
product of the one-electron Green's function G and the 
screened Coulomb interaction W)r- with a pseudopo- 
tential plane-wave basis. This GW/BSE approach has 
been implemented in the AI2NBSE and OCEAN pack- 
ages for valence and core-level spectra respectivelyJ^iii 
The method is advantageous compared to independent- 
electron approximations in that it provides a first- 
principles method for the inclusion of both quasi-particle 
and excitonic effects. These many-body effects modify 
peak positions, widths, and strengths, and hence they are 
crucial for a quantitative treatment and interpretation 
of the spectra. GW/BSE calculations have been carried 
out previously for both the valencei^ii^ and core^ spectra 
of water, with various approximations for the screened 



particle-hole interaction. A key difference is that our 
implementation uses a complex, energy- dependent GW 
self-energy, in addition to the screened particle-hole in- 
teractions. 

The remainder of this paper is as follows: The theoret- 
ical methods are summarized in Sec. II. Results for two 
forms of ice are presented in Sec. Ill, and for a model wa- 
ter system in Sec. IV. Finally, Sec. V contains a summary 
and suggestions for further work. 

II. THEORETICAL METHODS 

A. GW/BSE Approach 

The theory behind the GW/BSE approach imple- 
mented in AI2NBSE and OCEAN for valence and core 
excitations is described in detail in Refs. [13 and [ill re- 
spectively, so only a short summary is included here. The 
BSE represents the equation of motion of a particle-hole 
state, here an electron photo-excited into the conduction 
bands from either the occupied valence bands (optical) 
or a deep Oxygen Is orbital (x-ray). Briefly, the open- 
source plane-wave, pseudopotential code ABINIT— 
is used to calculate both occupied and unoccupied 
Kohn-Sham density functional theory (DFT) states of 
the ground-state Hamiltonian, which serve as a ba- 
sis for the BSE. This Kohn-Sham calculation uses the 
Ceperley-Alder, Perdew-Wang local density approxima- 
tion for the exchange-correlation potentialJi DFT or- 
bitals only strictly describe a non-interacting system, but 
they can be good approximations to quasi-particle wave- 
functions and perturbativley corrected, e.g., via Hedin's 
GW approximation. In this work, first-order quasi- 
particle energy corrections are added through the GW 
many-pole self-energy (MPSE) approximation of Kas et 
al.-- In contrast to previous GW/BSE approaches for 
these systemSfi^ii^ our MPSE is complex and energy- 
dependent, and thus accounts for both the energy shifts 
and inelastic losses which stretch and damp the spectral 
features. This MPSE model is based on a many-pole fit 
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to AI2NBSE calculations of the valence loss spectrum 
-lme-i(w). 

The optical and x-ray spectra, the latter including 
both x-ray absorption (XAS) and non-resonant inelastic 
x-ray scattering (NRIXS), which is also referred to as x- 
ray Raman spectra (XRS), are then calculated using the 
National Institute of Standards and Technology (NIST) 
BSE solver in AENBSEi" and OCEANii respectively. 
The core spectra in OCEAN use atomic core-level states 
and projector-augmented-wave (PAW) transition matrix 
elements--?- The BSE kernel includes both a screened di- 
rect attraction between the electron and hole and an un- 
screened repulsive exchange term. The treatment of the 
particle-hole interactions is an important consideration 
due to the weakly screened excitonic effects in water and 
ice. In both our valence and core codes the screened 
Coulomb interaction W is approximated as statically 
screened. A Hybertsen-Levine-Louie dielectric function 
is implemented in AI2NBSE, while OCEAN uses the 
random-phase approximation at short-range, switching 
to a model dielectric function at long rangei ^°i^^ In con- 
trast, previous approaches have used a variety of approx- 
imations ranging from a self-consistently screened DFT 
core-hole^'*^ to half-core-hole transition-state potentials)^ 



B. Structures 

The pair distribution function (PDF) g{R), and total 
electronic density of states (DOS) provide useful mea- 
sures of the ground-state structural and electronic prop- 
erties of our model ice and water systems (Fig. [T]). In 
order to interpret x-ray and optical spectra, it is im- 
portant that the model structures used in the calcula- 
tions match observed structural properties, e.g., the PDF 
from neutron or x-ray diffraction. For our model ice-Ih 
cell goo{R) is in reasonable agreement with experiment. 
However, the first shell peak is shifted to a slightly larger 
mean radius. To be consistent with experiment a root- 
mean-square disorder of about 0.25 A is needed, which is 
simulated in Fig. la by convolving the model PDF with 
a Gaussian. For our 17-molecule water cells (see Sec. Ill) 
fair agreement with experiment is observed, but with a 
slightly elongated first-shell 0-0 distance. The position 
and height of the first peak in goo{R) are sensitive to 
both the momentum range in the diffraction measure- 
ments and assumptions about the structure and core po- 
tentials. Discussions of these limitations can be found 
elsewhere, e.g., Refs.iGl and 

The total electronic density of states (DOS) provides 
a useful picture of the ground-state electronic structure, 
and in particular, the occupied and unoccupied energy 
levels in these materials. A comparison of the DOS 
for liquid, solid, and gas-phase H2O as calculated with 
ABINIT is shown in Fig.[l}). The DOS of both liquid wa- 
ter and ice-Ih are roughly similar, apart from broadening 
and the sharp peak at about -1-10 eV in ice. The similar- 
ity among the occupied DOS suggests that the core states 
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FIG. 1: (Color online) a) The oxygen-oxygen pair distribution 
function, goo{R), for our model structure unit cells (black, 
solid line) compared with experimental results for \ce^ (blue 
plusses) and two different water results (red boxes^^ and green 
crosses^^); b) comparison of the total electronic density of 
states for ice-Ih, liquid water and gas-phase models, calcu- 
lated using ABINIT. The broadening and vertical scaling are 
the same for all three structures, and energies are relative to 
the Fermi level. The features for the liquid cell are naturally 
broadened by structural disorder. 

are well localized and hence essentially molecular in char- 
acter. Likewise, the similarity in the unoccupied DOS 
between liquid and solid H2O suggests that short-range 
order dominates the spectra. The additional broadening 
in the liquid reflects the larger configurational disorder in 
amorphous structures. In contrast, the unoccupied DOS 
of the gas phase shows considerably more structure than 
either the liquid or solid phases. These results for the 
DOS are similar to previous work, e.g., Chen et ali^ 



C. GW Self-energy 

Quasi-particle self-energy effects are of crucial impor- 
tance in broad-spectrum calculations of optical and x-ray 
spectra, because independent-particle calculations sys- 
tematically underestimate conduction band widths. A 
comparison of our GW many-pole self-energies for wa- 
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FIG. 2: (Color online) Comparison of GW self-energies for 
solid and liquid water as calculated with our MPSE model. 
Note that these are quite similar, indicative of similar broad- 
ening and self-energy damping effects in water and ice. The 
quasi-particle self-energy broadening may be compared to 
other sources of broadening in the spectra: experimental reso- 
lutions (« 0.5 eV), the O Is core- hole lifetime (0.16 eV), and 
vibrational broadening which is yet smaller and not included. 



ter and ice is shown in Fig. (3] Note that the imaginary 
parts are quite similar and grow significantly above about 
15 eV, while the real parts differ little. This reflects the 
close similarity between the valence dielectric response 
(i.e., £2 and — Ime^^) of water and ice, and leads to an 
overall stretch of the spectra by about 1 eV. Unlike full 
GW approximations however, the MPSE fails to fully 
correct the gap in these systems, and the optical spectra 
that follow have been aligned to match the observed gap. 



III. ICE SPECTRA 

Water has a rich phase diagram with at least 12 crys- 
talline phases of ice, 9 of which are stable, in addition 
to several amorphous phases. Complicating the discus- 
sion of H2O, ice is proton disordered, but this cannot 
be treated exactly by our approach which enforces peri- 
odic boundary conditions. Of the known, stable phases 
of ice, only ice- VIII has proton ordering with an anti- 
ferroelectric unit cell. Recent investigations of the pro- 
ton ordering phase transition between ices VII and VIII 
found no significant difference in the Oxygen K-edge 
NRIXS, implying that the effect of proton disorder on 
that spectrum is small. 

In addition to local fluctuations in the dipole moments 
of most phases of ice, the low-mass hydrogen atoms have 
a significant amount of zero-point motion on top of ther- 
mal vibrations, which is not captured by standard Born- 
Oppenheimer molecular dynamics. Once again, accurate 
treatment of this disorder requires large unit cells which 
are not feasible in our approach. Recent work has focused 
on capturing the quantum nature of the hydrogen move- 
ment, but such improvements are not included here 3^ 





AI2NBSE + MPSE 




'. AlilNDot 




I Expt. 


















J 

















5 10 15 20 25 30 

Energy (eV) 



FIG. 3: (Color online) The imaginary part of the dielectric 
constant of ice-lh calculated using AI2NBSE with the MPSE 
correction (solid red), and with a static energy shift (dashed 
green). For comparison an experimental measurement of ice- 
lh is shown (dotted blue).— The discrepancy likely reflects a 
lack of disorder in our model structure (see text). 



A. Ice Models 

For this study we have focused on two forms of ice; 
ice- VIII, which is stable below 273 K and pressures in 
the range of 2 GPa to 50 GPa, and ice-lh, which is the 
common form of ice, stable under ambient pressure. Our 
ice-lh model uses a 16-molecule cell determined by en- 
forcing the experimental density and following the ice 
rules for hydrogen placement with an 0-H bond distance 
and angle of 1.01 A and 106° respectively. A realistic 
simulation of ice-lh should have full proton disorder, but 
our models are based on small cells with periodic bound- 
ary conditions, and thus contain artificial order. Further 
calculations for other ice-lh models using larger unit cells 
and finite temperature configuration sampling are called 
for to assess the validity of our simplified model. 

Ice- VIII is formed under high pressure at temperatures 
below the conventional freezing point of water. The pres- 
sure requirements preclude experimental techniques in- 
compatible with diamond-anvil cells such as lower energy 
XAS measurements. However, NRIXS can probe the 
same final-state as XAS with an additional tunable pa- 
rameter q, the momentum transfer and at a much higher 
beam energy. The structure of ice- VIII is that of two in- 
terpenetrating cubic ice lattices resulting in a compressed 
second shell and higher densities than ice-lh. The unit 
cell contains only 4 molecules and has anti-ferroelectric 
Hydrogen ordering, making it an easier system to simu- 
late than standard ice. Our ice- VIII cell uses the experi- 
mental lattice constants at a pressure of 2.4 GPai^ 



B. Ice valence spectra 

For the ice-lh spectra, the wavefunctions from ABINIT 
were calculated on a 4 x 4 x 4 fc-point mesh with 464 bands 
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FIG. 4: The calculated momentum-transfer dependence of 
the oxygen K-edge in ice-Ih (solid red), including the MPSE. 
For comparison experimental results (dashed blue) for XAS"^^ 
and NRIXS2^ are shown. 

(400 conduction bands). For comparison, the struc- 
ture and valence spectrum of hexagonal ice (ice-Ih) were 
studied previously with DFT/BSE approaches, yielding 
qualitative agreement with experiment.— For the valence 
spectrum (Fig. , the calculated excitonic peaks at 8 eV 
and 10 eV are significantly stronger than those measured. 
Between 10 eV and 20 eV the absorption is also too 
strong and too thinly peaked, compared to experiment 
which exhibits additional broadening past 15 eV. These 
discrepancies are partly due to the lack of disorder in our 
model structure, and point to the need to consider more 
elaborate models that can account for such disorder. The 
effect of the MPSE is seen as a stretch of the spectra, to- 
gether with energy dependent broadening. Coupling be- 
tween electronic excitations and phonons is expected to 
provide small additional broadening, but is not included 
in this work. 



C. Ice core spectra 

For the oxygen K-edge XAS and NRIXS calculations a 
4x4x4 fc-point mesh and 400 conduction bands were used 
for the BSE calculation and 900 bands for the screen- 
ing. Oxygen K-edge XAS and NRIXS results for ice-Ih 
are shown in Fig. |4l All calculated x-ray spectra have 
been broadened by convoluting with a Lorentzian to ac- 
count for the core-hole lifetime and a Gaussian to match 
reported experimental broadening. Further damping is 
provided by the imaginary part of the self-energy (cf. 
Fig. [5]). Results are shown averaged over orthogonal in- 
cident photon polarizations or momentum transfers for 
XAS and NRIXS respectively, since experimental data is 
reported for polycrystalline samples. 

Our XAS calculation matches experiment fairly well, 
but has a lower pre-edge (535 eV) and is noticeably too 
narrow in overall width (Fig. U]) . Recently the momentum 
dependence of the O K-edge in ice has been measured us- 
ing NRIXS.25 The theoretical simulation of NRIXS with 
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FIG. 5: (Color online) The calculated (solid red) momentum 
transfer dependence of the oxygen K-edge in ice- VIII NRIXS 
compared with experiment (dashed blue)'^^ for small {q = 
4.2 A~^) and large {q = 9.4 A~^) momentum transfer. 



OCEAN differs only slightly from XAS in that the fi- 
nite momentum transfer q can break the dipole selection 
rules, allowing transitions to s-type and higher angular 
momentum final states. Overall the calculations yield 
fair agreement for the spectra, including both the rela- 
tive weights and the momentum dependence of the pre- 
edge feature (Fig. |4]). However, the balance between the 
s- and p-character of the pre-edge is shifted too much 
towards the s-type in our calculation, leading to overly 
strong growth of the pre-edge with increasing q. Note 
that experiment shows little difference between g « 
(XAS) and q = 3.1. 

We have also examined the momentum dependence of 
ice- VIII (Fig. [5|) which exhibits a similar evolution with 
increasing q as ice-Ih, and good general agreement with 
experiment. Like ice-Ih, our calculations of the O K-edge 
in ice- VIII are too strongly peaked in the main edge and 
have too little weight above about 543 eV when com- 
pared to experiment, but show improvements over pre- 
vious theoretical results.— Similar to ice-Ih, the calcu- 
lations of ice- VIII also show excessive q dependence for 
the pre-edge, thus suggesting either a limitation of the 
DFT wave functions or the need to consider thermal and 
quantum disorder. 



IV. WATER: 17-MOLECULE LIQUID CELLS 

A. Water structural model 

Our calculations for optical and x-ray spectra were 
carried out on 17-molecule snapshots for liquid water 
obtained from the molecular dynamics (MD) results of 
Garbuio et alJ^ Despite the small size of these model 
structures, they already exhibit substantial disorder and 
are found to have a pair-distribution function (Fig. [1]) in 
fairly good agreement with experiment. However, the O- 
O nearest neighbor distance is notably too long by a few 
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FIG. 6: (Color online) Calculated 62 spectra for a single MD 
snapshot with a 17-molecule water cell using the MPSE (solid 
red) and a static energy shift (dashed green) compared to 
experiment (dotted blue).— 

tenths of an angstrom. Additionally, while bulk water 
is a disordered system with no net dipole moment, the 
limited cell size of our systems precludes this constraint; 
for our samples the moments were found to vary widely 
from 14 Debye to 28 Debye. 



B. Water valence spectra 

For the valence calculation a 4 x 4 x 4 /c-point mesh and 
300 conduction bands were used, and both a static energy 
shift of 4.13 eV, as suggested by Garbuio et al.,!^ as well 
as the energy dependent MPSE were applied. With the 
static shift the positions of the initial peaks match well 
with experiment (Fig. ^ and in general, the agreement 
is excellent over the full range of the spectra. Our results 
for the peak positions in €2 agree qualitatively with those 
of Garbuio et al., but yield features that are smaller in 
overall magnitude. The energy dependence of our MPSE 
(Fig. [2]) introduces a stretch in the spectrum of w 5 % in 
the range of 10 eV to 20 eV, while Garbuio et al. reported 
the GW energy correction to be almost constant across 
the low-lying conduction bands for water. 



C. Water core spectra 

XAS calculations were carried out for 8 of the MD 
snapshots using a 4 x 4 x 4 fc-point mesh including 300 
bands for the final states and 400 conduction bands 
for the screening. The screened core-hole potential was 
found to be nearly identical for different O sites within 
a cell, reflecting the molecular character of liquid water. 
The O K-shell XAS calculation (Fig. O exhibits consid- 
erable variation among the different oxygen sites both 
within a cell and between cells. Both the calculated and 
recent experimental spectra'^^''^"^ exhibit a notable inten- 
sity shift and increased broadening going from ice to wa- 



FIG. 7: (Color online) The XANES spectra for the O K- 
edge in water. The OCEAN results, including the many-pole 
self-energy, are compared with recent Scanning Transmis- 
sion X-ray Microscopy (STXM>22, and Total Electron Yield 
(TEY) experiments.-^ The error-bars for the calculation rep- 
resent the rms variation between the 8 different MD snap- 
shots, showing reasonable similarity between different geome- 
tries sampled by the MD despite the limited cell size. Inset; 
The main edge for the MD snapshots with the lowest and 
highest dipole moments compared to average (solid). With 
decreasing dipole moment the first peak intensity decreases 
with a slight shift to higher energies, a trend consistent with a 
decrease in the discrepancy with experiment as the net dipole 
moment approaches zero. 

ter, shifting spectral weight to the lower part of the spec- 
trum from 534 eV to 538 eV. The calculated O K-edge 
spectra is in good agreement with experiments when cor- 
rected with the MPSE, but, as for the calculated ice-Ih 
spectra, the main edge at about 536 eV is stronger than 
in experiment. The peak at 541 eV in experiment is ev- 
ident in the calculation but is lacking in strength. The 
success of the MPSE in correcting the overall width of the 
near-edge is consistent with calculations using the COH- 
SEX approximation,'* but in contrast to earlier work— 
which showed that the self-energy for water was nearly 
constant for their G^W'^ approximation. 

In a disordered system like liquid water, each oxygen 
is subject to a different static potential and to slightly 
different core-hole screening by the valence electrons, re- 
sulting in shifts in the binding energy of the Is electrons. 
Exact, absolute energies were not calculated, but relative 
shifts were determined according to the relation 

Eu - E'^: + y^''^ + ^Wc. (1) 

-E™ is the site-independent binding energy which has 
been aligned to match experiment. The second term 
V^^ , is the total Kohn-Sham potential at the site, in- 
cluding the effects of all the other ionic cores. The last 
term is the effect of the spectator electrons screening the 
core hole, reducing the energy necessary to excite the 
Is electron. Both V^^ and Wc are evaluated at the 
Oxygen site. These binding-energy shifts lead to a slight 
broadening of the spectra in Fig. [7] 
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We also find that the unphysical dipole moments of 
the small unit-cells in this study have a systematic effect 
on the calculated spectra (inset of Fig. [7]). That is, the 
two main peaks of the calculated O K-edge at 537 eV and 
540 eV exhibit shifts in intensity and position that loosely 
correlate with dipole moment, i.e.,. the local, static elec- 
tric field felt by each water monomer in the ground state. 
At higher values of the dipole moment the main edge fea- 
tures shift to lower energies, increasing discrepancy with 
experiment. This again points to the need for simula- 
tions with larger cells that more accurately characterize 
the model structure and properties. 

V. SUMMARY AND FUTURE PROSPECTIVE 

Theoretical calculations have been carried out for both 
valence and core x-ray spectra of a number of model ice 
and water systems using modern GW/BSE theoretical 
methods. Our results suggest that accurate calculations 
of the spectra for these systems are quite sensitive both to 
the theoretical methods and details of the model struc- 
tures. Nevertheless, we have shown that the improved 
treatment of many-body effects in the present approach 
yields better agreement with experiment in terms of both 
relative peak weights as well as overall width and feature 
locations compared to calculations that ignore these ef- 
fects. 

We find that the inclusion of an accurate quasi-particle 
self-energy is important to characterize the damping and 
self-energy shifts in the spectra. In particular we find 
that the stretch provided by our GW MPSE significantly 
improves the agreement between calculations and exper- 
iment for the ice and water systems. Our results for the 
valence £2 spectra for the 17- molecule water model dif- 
fer somewhat from those of Garbuio et al. Theirs have 
somewhat larger values of £2, while ours exhibits an ex- 



tra excitonic peak. The origin of these differences is likely 
due to differences in the screening approximations used. 
In any case, our present results appear to be in reason- 
able agreement with recent experimental results both for 
the core- and valence-spectra of water. However, we sug- 
gest that configurational averages should be carried out 
for ice-Ih and ice- VIII to better understand the effects 
of disorder, finite temperature, and zero-point variations 
in structure. A greater variation in hydrogen position- 
ing and bonding from finite temperature and zero-point 
motion effects could lead to changes in feature weights 
as well as adding an increased broadening from disor- 
der. Additionally, convergence should be checked with 
respect to the size of the cells being used, especially be- 
cause larger cells should allow for better control of the net 
dipole moment. At present, however, system size is lim- 
ited by our codes and hence significantly larger cells are 
computationally impractical. The oxygen K-edge should 
be calculated for a sufficient number of MD snapshots to 
ensure a good representation of the physical system. Site- 
specific geometrical information, i.e., bond length, an- 
gles, and hydrogen bond coordination, can then be com- 
pared to contributions from individual oxygen XANES to 
correlate differences in local environment and calculated 
XAS/NRIXS. 
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